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I. INTRODUCTION 

The unintegrated parton distributions (UPD's) have 
been considered in numerous works on applications of 
perturbative QCD QSSSIEIESIEIa E3, EH E3, 
El llil [TEl ITi| . These distributions are generalizations 
of the usual integrated parton distributions (PD), and 
in some sense more basic objects, as PD's are obtained 
from UPD's by integration over the transverse momen- 
tum of the parton. The notion of UPD's relies on the k±- 
factorization, and, in the spirit of the CCFM equations 
[ttI ITsl . flfl l20j . introduces two scales: the probing scale, 
Q, and the transverse-momentum scale, k±. Recently, 
the UPD gained substantial attention, since they enter 
many exclusive physical processes, such as t he p roduc- 
tion of the Higgs boson B^, > the W boson |23 | . heav y 
flavors HI H |H E^Jpa, the jet production[2ll Hf, 
particle production 29], or hadron production in rela- 
tivistic heavy-ion collisions [30l l3l| . The unintegrated 
distributions were also used in studies of the longitudinal 
[32| , charmed j3^| , and spin 0, [35| structure functions 
of the nucleon, as well as analyzed in the dipole picture 
of QCD [13. 

The UPD's, similarly to other entities in QCD, un- 
dergo evolution with the change of the probing scale Q. 



'Dedicated to the memory of Jan Kwiecihski 

T Electronic address: earriola@ugr.es 

t Electronic address: Wojciech.Broniowski@ifj.edu.pl 



The philosophy adopted here is similar to the case of the 
integrated PD. At some initial scale Qo we need to know 
the non-perturbative quantity from measurements, mod- 
els, or lattice calculations, and then we can evolve it to 
a different scale Q with the help of suitable QCD evolu- 
tion equations. In the case of integrated PD we need to 
assume the dependence of PD on the Bjorken x variable 
at the initial scale Qo. For the UPD we need to know 
in addition the dependence on the transverse momentum, 
k±, or, equivalently, the transverse coordinate b, which is 
the variable Fourier-conjugated to fcj_. Knowing this, we 
compute, with no extra physical input apart for the as- 
sumptions entering the QCD evolution, the unintegrated 
distribution at the final scale Q. 

Important physics questions may be answered. In par- 
ticular, the UPD evolve in such a way that the average 
transverse-momentum increases in a specified way with 
the scale [lj; UM Gil HH • This spreading can be stud- 
ied quantitatively within the approach. This constrains 
the freedom in phenomenological analyses of processes 
involving the UPD's. 

In his studies of the problem . Kwiecihski started from 
the CCFM formalism [13 El E HI, which explic- 
itly involves two separate scales: the probing scale, Q, 
and the transverse momentum of the parton, k±. The 
CCFM equations were subsequently extended to include 
the quarks, as well as reduced to the single- loop approx- 
imation. In addition, the non-Sudakov form factor was 
dropped. Kwiecihski found that in this approximation 
the equations diagonalize in the space Fourier-conjugated 
to the transverse momentum, where they assume a par- 
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ticularly simple and eleg ant form in a close resemblance 
of the DGLAP [H HalM |H equations for the inte- 
grated PD. The only, but most important, difference is 
the appearance of the Bessel function Jo(Qb) in the evo- 
lution kernel. Thus, the evolution depends on the trans- 
verse coordinate b. In the original work of Ref. ^| these 
equations were called "the CCFM equations for the UPD 
in the transverse-coordinate space in the single-loop ap- 
proximation" . Due to numerous steps leading away from 
the original CCFM, we find it more appropriate to call 
the equations of Refs. [H| the Kwieciriski equations for 
the UPD evolution. Since in the case of integrated PD 
these equations reduce to the usual LO DGLAP, the 
range of applicability of the Kwiecihski equations is not 
larger as for the LO DGLAP, with not too small and not 
too large values of the Bjorken x variable. 

Formally, the Kwiecihski equations are integro- 
differential equations with the kernel depending on the 
transverse coordinate b (cf. Sec. [HJ. As such, they 
are not trivial to solve numerically. The method of the 
original works [iH ITHl ITfjj involved the Chebyshev inter- 
polation in the x and Q spaces for each value of b. In this 
paper we offer an alternative method, based on the evo- 
lution of the x-moments of UPD. In the moment (Mellin) 
space, the evolution equations become a set of ordinary 
differential equations, which can be solved numerically 
in a very efficient way (see Sec. Illlf) . We derive analytic 
expressions for the 6-dependent anomalous dimensions, 
which can be written in terms of hypergeometric func- 
tions. Then, the inverse Mellin transform to the original 
x-space is performed via numerical integration with os- 
cillatory functions. We show that the procedure is fast 
and stable, providing a useful numerical tool for evolving 
the UPD. 

The Mellin-transform method allows us to carry ana- 
lytic considerations, such as studies of certain limits of 
the equations, specifically the cases of low and and large 
b, and x approaching the endpoints. We pursue these 
considerations, which can be done since the form of the 
6-dependent anomalous dimensions is analytic. 

In addition to developing a different numerical method 
fSec. lIIIf) . our study differs from Ref. jig in two physical 
aspects. Firstly, rather than guessing the initial shape in 
b, we use the results of low-energy chiral quark models 
(Sec. II V|) . We consider two models: the recently pro- 
posed Spectral Quark Model of Ref. j illlilli^ and the 
Nambu-Jona-Lasinio model with the Pauli Villars regu- 
larization . These models were used before to describe 
the integrated PD EH 113 El] , and were shown to do 
a surprisingly good job, in particular for the valence dis- 
tribution in the pion. They were also used to describe 
successfully other aspect of high-energy processes, such 
as the pion distribution amplitude [491 and generalized 
(off- forward) PD of the pion [5(J, |51| • The models give 
the initial condition at the model scale Qq in a partic- 
ularly simple, factorized form. The valence quarks are 
distributed uniformly in x, while the gluons and the sea 
quarks vanish. The 6-dependence is a simple, analytic 



function with exponential fall-off at large b. We stress 
that the b dependence is a prediction of the model, rather 
than a mere guess, as is frequently made. Secondly, our 
implementation of the evolution switches from three to 
four flavors above the charm-production threshold, cus- 
tomarily taken at Q 2 — 4 GeV 2 . Our numerical results 
are presented in Sec. where we show the dependence 
of the UPD's on x and b. 

Since the analytic forms of the anomalous dimensions 
involve generalized hypergeometric functions, which may 
be cumbersome to program, we have developed a low-& 
expansion (Sec. IVIfl . which makes the calculations sim- 
pler when b is not too large. The expansion is in powers of 
^Aqcd, and is fast and stable. For the opposite problem, 
where Qb is large, we have obtained asymptotic expan- 
sions (Sec. rVTQ . which allow to deal numerically with the 
generalized hypergeometric functions. 

Our method of solving the equations in the Mellin 
space carries additional bonuses. In particular, it allows 
us for analytic considerations in the investigation of for- 
mal limits at x — ► 1 (Sec. IVlllj) and x — > (Sec. lDT)l . 
At x — > 1 we show that the 6-dependent non-singlet dis- 
tribution approaches very fast the integrated non-singlet 
distribution. At x — > we find generalizations of the 
double-leading- logarithm (DLLA) formulas of Ref. [54| . 
Finally, we examine the large-6 behavior, where we show 
that the evolution-generated UPD's from the Kwiecihski 
equations exhibit power-law behavior at large b. The 
fall-off is much faster for the gluons than for the quarks. 

Widening in the transverse momentum of all partonic 
distributions is confirmed. We show that (k\) grows 
with the probing scale as Q 2 a(Q 2 ). We write an approx- 
imating formula for the width for each partonic species, 
which may be useful in practical applications with the 
pion (SeaEJ. The widening effect becomes stronger and 
stronger as Q increases or x decreases, and it is bigger 
for the gluons than for the non-singlet and singlet quarks 
(see Sec. Ixl) . 

The numerical method of this paper, which is easy to 
program and numerically fast and efficient, can be used 
for other initial conditions as well, for instance for the 
GRS parameterization of the pion or the GRV pa- 
rameterization [5^], of the nucleon, supplied with a pro- 
file in the transverse coordinate, as originally studied in 
Ref. ^(|- The only difference is in the form of the ini- 
tial Mellin moments, which acquire the dependence on 
b. General prediction of the method in formal limits are 
listed in Sec. 

Appendices contain many technical details, such as 
the perturbative QCD parameters and splitting functions 
(App. 0, the analytic form of the 6-dependent anoma- 
lous dimensions which enter the evolution in the Mellin 
space (App. [Bj, their low-6 (App. [UJ and high-& ex- 
pansion (App. [Dj, as well as the pole-residue expansion 
(App. [EJ. The latter is useful in analytic considerations 
near x = 0. 



3 



II. THE KWIECINSKI EQUATIONS 

In his studies of the UPD's, Kwiecihski El El EH El 
started from the CCFM formalism 

El El El la explic- 
itly involving two separate scales: the probing scale, Q, 
and the transverse momentum of the parton, k± . Then, 
the original CCFM equations were supplemented with 
the quarks, as well as reduced to the single- loop approx- 
imation. The latter approximation replaces the angular 
ordering of the emitted gluons with the ordering of their 
transverse momenta. In addition, the non-Sudakov form 
factor was dropped. Kwiecihski realized that the evolu- 
tion equations for the UPD's acquire a particularly sim- 
ple form in the transverse-coordinate space, 6, conjugated 
to the transverse momentum k±. For each distribution 
one introduces 

/>oo 

fj{x,b,Q) = / tordk±k± J (&fcx) fj{x, k±,Q), (2.1) 

where j = NS (non-singlet quarks), S (singlet quarks), 
or G (gluons), and Jo is the Bessel function. In order to 
avoid confusion, we stress that the transverse coordinate 



b, conjugated to the parton's transverse momentum, is 
not the impact parameter, appearing in the analysis of 
the generalized PD. The latter quantity is conjugated the 
the transverse momentum transfer in off-forward scatter- 
ing processes. 

At b = the functions fj are related to the integrated 
parton distributions, pj(x,Q), as follows: 

f j (x,0,Q) = ^Pj(x,Q). (2.2) 

More explicitly, for the case of the pion studied in this 
paper (we take n + for definiteness) we have 

Pns — u — u + d — d, 

p s = u + u + d + d + s + s + (2.3) 

Psca = ps -PNS = 2d + 2u + s + s + 

PG = 9, 

where . . . stand for higher flavors. 
The Kwiecihski equations read |l6j : 



Q 



2 dfas{x,b 1 Q) 
dQ 2 
2 df s (x 7 b,Q) 



dQ 2 



Q 



2 df G (x,b,Q) 



dQ 2 



dzP qq (z) 



e 



a s (Q 2 ) 
2tt 

J dz(e(z - x) J ((l - z)Qb) 

[zP qq (z)+zP Gq {z)]f s {x,b,Q) 
« S (Q 2 ) '•' 



2tt 



(z - x) J ((l - z)Qb) f NS b, Q) - f m (x, 6, Q) 

P qq (z)fs{^,b,Q)+P qG (z)f G (^,b,Q) 



P Gq (z) fs (-, 6, 0) + Pgg(z) fo (-, 6, Q) 



dz<e(z-x)J ((l-z)Qb) 



[zP GG {z) + zP qG (z)]f G (x,b,Q) 



r 



(2.4) 



The splitting functions P a b{z) are listed in Eq. (|A9(I . 

Following Ref. E3j a factorized form of the distribu- 
tion functions at the initial scale Qq is assumed, 



F 



NP 



{ty^PiiX'Qo), 



(2.5) 



with the profile function F NP (b) taken to be universal 
for all species of partons. The factorization assumption 
(|2.5|) is technical and one can easily depart from this 
limitation in numerical studies. We note, however, that 
the models of Ref. 

sum, 

studied in Sec. II VI do predict 
a factorized initial condition of the form (|2.5|) . The input 
profile function, F NP (6), is linked through the Fourier- 
Bessel transform to the k± distribution at the scale Qo. 
At b = the normalization is F NP (0) = 1. The profile 
function factorizes from the evolution equations. This is 
clear, since any solution of Eq. I|2.4|l remains a solution 



when multiplied by an arbitrary function of 6. Due to 
evolution, at higher scales Q we have 



fj(x,b, Q) 



F* p (b)f? m \x,b,Q), 



(2.6) 



with fj Vol (x,b,Q) satisfying equations identical to (12.41) . 

We should stress again an important physical differ- 
ence between F NP (6) and /f ol (x, b, Q). While F NP (b) 
originates entirely from low-energy, non-perturbative 
physics, fj Vol (x, b, Q) is given by the perturbative QCD 
evolution with Eq. (|2.4|l from the initial condition 



ff°\xAQo) = -^Pj(x,Qo). 



(2.7) 



Throughout this paper, except for Sec. IIVI we fo- 
cus on the perturbative evolution and the functions 
/ cvol (x,6, Q). The function F NP (b) is referred to as the 
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initial profile and /j (x, b, Q) as the evolution- generated 
UPD. 



In the singlet channel we find the coupled set of equa- 
tions, 



III. THE KWIECINSKI EQUATION IN THE 
MELLIN SPACE 

We define the x-moments of an unintcgratcd parton 
distribution in the impact parameter space as (we retain 
the same symbol for the function and its Mellin trans- 
form, hoping the distinction made by the argument pre- 
vents any confusion) 



fj(n,b,Q) = / dxx^f^^Q). 



(3.1) 



In the Mellin space, the evolution equations for the UPD 
are very simple, as they become diagonal both in b and n. 
They involve ^-dependent anomalous dimensions, equal 
to 

ln,ab(Q b ) = 4 / dz [z n J ((1 - Z )Qb) - 1] P ab {z) 

Jo 
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Inlb -4 / dz z n [Jo ((1 - z)Qb) - 1] P ab (z), 
Jo 



(3.2) 



where the values at b = are 

<7 (0) 



~(°) 
in,qq 



-Y (0) 
<n,qG 



-Y (0) 
ln,Gq 



!n,GG 



-4 / dz{z n ~l)P qq {z), 
Jo 

-4 [ dz [(z n - z)P qq - zP G , 
Jo 

-4 f dzz n P qG , 
Jo 



(3.3) 



-4 / dzz n P Gq , 
lo 



dz [{z n - z)P ( 



GG 



zP, 



qG\ 



Their explicit form for various channels is listed in 
Eq. l|B6IB7j) . The fact that we can write the analytic 
form of the integrals in Eq. (|3.2[) (see App. [BJ) allows 
for efficient numerical calculations and analytic consider- 
ations. 

The integration of both sides of Eq. I|2.4[) with 
j} dxx n_1 yields for the non-singlet case the equation 

dhs(n,b : Q) a S {Q 2 ) , , n , ,„ 
dQ 2 = — g n Q2 ln,Ns{Qb)fNS(n, b, Q). (3.4) 

The formal solution of Eq. H3.4|l can be readily obtained 
as 



hs(n,b,Q) 
hs(n,b, Q ) 



exp 



Q 
Ql 



dQ' 2 a{Q' 2 ) . , n . 
7i — 7ns [n, b, Q ) 



df s (n,b,Q) _ a S (Q 2 ) 
dQ 2 8ttQ 2 
x hn,qq(Qb)fs(n, b, Q) + ~fn, q G{Qb)f G (n, b, Q)] , 

df G {nAQ) _ as(Q 2 ) 
dQ 2 8nQ 2 
x hn,Gq(Qb)fs(n, b, Q) + j n ,GG(Qb)f G (n, b, Q)] , 

(3.6) 



which has the formal solution 



fs(n,b,Q) 
fo(n,b,Q) 



(3.7) 



Vexp 



T n (Qb) 



Q2 dQ' 2 a(Q' [ 
x SttQ' 2 



-r»(Q6) 



fs(n,b, Q ) 
f G (n, b, Q ) 



ln m (Qb) ln,qG(Qb) 
ln,Gq{Qb) ln,GG{Qb) 



The symbol V indicates that powers of T n are ordered 
along the integration path. The qualitative difference 
between Eq. I|3.7|) and the LO DGLAP equations is the 
fact that r„ depends on the evolution variable Q. This 
makes the singlet sector more difficult to analyze analyt- 
ically. Equations (|3.(il) can be solved numerically for any 
value of n real or complex [5j| . For the case of integrated 
PD (b = 0) Eq. I|3.4|l and 13. 6|) reduce to the well-known 
LO DGLAP equation in the Mellin space. 

The corresponding UPD in the x space can be recon- 
structed using the inverse Mellin transform, 



fj(x,b,Q) 



no+ioc 



dn 



x- n f 3 {n,b,Q), (3.8) 



where uq has to be chosen in such a way as to leave all 
the singularities on the left-hand-side of the contour. It 
turns out that for b ^ the analytic structure of the 
b— dependent anomalous dimension remains the same as 
for the 6 = case. This can be inferred directly by 
studying the analytic structure of the formulas (|B6IB7|) 
or the pole-residue expansion of Eq. {E4j. Thus, we have 
the result that that "f n ,Ns{Qb), ln,qq(Qb), and ~f n , q G{Qb) 
have poles at n = —1,-2,-3,..., while j n ,GG(Qb) and 
ln,Gq{Qb) have poles at n = 0,-1,-2, .... Parameteriz- 
ing the contour in Eq. (|3.8|l as n = no + it, we arrive at 
the inversion formula 

f°°dt 

fj(x, b, Q) = x- no / — (cos(Hogx)Re[/, (n + it, b, Q)] 
Jo ^ 

+ sin(t\ogx)lm[fj{n +it,b,Q)\). (3.9) 

For the non-singlet case we take no = 0, while for the 
singlet case n = 1. As an additional check we have also 
verified that a bended integration path n = 



(3.5) with < r < oo, as used in [55|. also works nicely. 



re 



iir/4 
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IV. THE INITIAL CONDITION FOR THE UPD 
OF THE PION 

As an illustration of our method, we consider the 
UPD's of the pion with the initial condition at Qo pro- 
vided by two large-iV c low-energy chiral quark models. 
The original work of Ref. ^(| tested the GRS parame- 
terization for the pion |52j and the GRV parameteriza- 
tion for the nucleon [5^ , supplied with a Gaussian profile 
F NP (6). The considered models generate the functions 
F NP (b) as genuine model predictions of low-energy non- 
perturbative physics, with no freedom involved. Our first 
model is the recently proposed Spectral Quark Model 
fSQMIpll IH. |43| . and the second one is the popular 
Nambu-Jona Lasinio model (NJL) with the Pauli-Villars 
regularization |44| , treated for simplicity in the strict chi- 
ral limit. 

Firstly, since the chiral quark models have no gluon 
degrees of freedom, we have at the model scale 



or 



g(x,Q ) = 0, 



rj° i (x,b,Q a ) = o. 



(4.1) 



(4.2) 



For the integrated valence quark PD both models predict, 
in the chiral limit and at the model scale Qo, that 



q(x,Q o )=0(x)O(l-x) 



(4.3) 



i.e. a constant value with the proper normalization 
and support. Thus the corresponding J ovo1 functions of 
Eq. I|2.2[l are linear in x, 



f c N %(x,b,Q ) = x9(x)9(l-x), 



(4.4) 



The scale of the model, as found from the momentum 
sum-rule, is rather low: Qo = 313 MeV. Although this 
is admittedly a very low scale, one may hope that the 
typical expansion parameter, a(Qg)/(27r) ~ 0.3, is low 
enough to make the perturbation theory sensible. This 
claim gains support from the next-too-leading analysis of 
the integrated PD jl^, where the corrections are found 
to be small. 

The QCD evolution is crucial for the phenomenological 
success of the considered low-energy chiral quark models. 
In Refs. [ill l46l l47j it has been found that the non- 
singlet distribution, when evolved to the scale of 2 GeV, 
agrees very well with the SMRS parameterization of the 
pion data |56j . while in |57j it has been very favorably 
compared to the old E615 data at 4 GeV [5£| (see the 
discussion in Sec. El and Fig. [5j|. 

In the SQM, the valence UPD of the pion at the scale 
Qo is 



q(x, k±,Q ) = q(l-x,k±) 



(4.5) 



6M£ 



ir(k 2 ± + M 2 /4) 5 / 2 



9(x)9(\ - a;), 



where My = 770 MeV is the mass of the p meson. 
Passing to the impact-parameter space with the Fourier- 
Bessel transform yields 



q(x,b,Q ) = 2n k±dk±q(x, k±)J (kj_b) (4.6) 
Jo 



bA'h 



exp 



The expansion at small b gives 



q(x,b, Qo) 



1 



mU 2 



M v b 



Myb 3 



6(x)6(l - x). 



8 24 
x9(x)9(l -x), 



(4.7) 



and average transverse momentum squared is equal to 

4 dq(x, b) 



(k 2 nSQM _ Jd 2 k_L k 2 ± g(x,k±) = 

[ ±;NP ~ Jd 2 k ± q(x,k ± ) q(x,b) db 2 



M 2 



b=0 

(4.8) 



which numerically gives (fci)j^ M = (544 MeV) 2 (all at 
the model working scale Qo). The subscript NP re- 
minds us that the quantity comes entirely from the non- 
perturbative physics, entering profile function i 7 ' NP (6) 
(see the discussion at the end of Sec. . 

In the NJL model with the PV regularization the 
analogous formulas read |4Sj 



q(x, k±, Qo) = q(l - x, k±) 
A 4 M 2 N r 



4/ 2 7r 3 {k]_ + M 2 ) {k]_ + A 2 



M N / 

q(x,b,Q ) = —. 2-f (2K (bM) ~ 2Ko(bVA 2 +M 2 ) 



M 



n 



(4.9) 
-9(x)e(l-x), 



4/ 2 

bh 2 Ki (by/ A 2 + M 2 ) * 
VA 2 + M 2 



(x)6(l-x), 



(4.10) 



where the pion decay constant is given by 

M 2 N C ( A 2 + (A 2 + M 2 ) log Af±^ 

f 2 = ^ ' . (4.11) 

U 4tt 2 (A 2 + M 2 ) V ' 

The parameters of the model are adjusted in such a way 
that U = 93 MeV, namely M = 280 MeV and A = 
871 MeV. The expansion at small b yields 

q(x,b,Q )= (4.12) 
' Af 2 V c (A 2 -Af 2 log^±^)& 2 

1 le^T 2 

x9(x)9(l-x), 
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FIG. 1: The initial profile functions Fnp for the SQM and 
NJL models, plotted as functions of the transverse coordinate 
b. The fall-off is exponential, according to Eq. (14. 14|l and 
KWt . 



and average transverse momentum squared is equal to 
M 2 N c {k 2 -M 2 logA!±*£l 



2 \NJL 
NP 



^ 2 f 2 



(4.13) 



which numerically gives {kj_)^p L — (626 MeV) 2 (at the 
scale Qo)i which is similar to the number from the SQM. 

Finally, in the notation of Eq. i|2.6[) we can write that 
the initial profile function is 



R 



sqm(^) — 
njlO) 



= 1 



bMy 



exp 



Mvb 



(4.14) 



M 2 N. 



-| (2K {bM) - 2K a (bVA 2 +M 2 ) 

bA 2 K 1 (bVA 2 + M 2 ) \ 

VA 2 + M 2 )' [ ' ' 

Both initial profile functions are displayed in Fig.Q] Note 
that the profiles, although have a b 2 correction at small 
b, are not Gaussian, and at large b display an exponential 
fall-off. 

As discussed in Sec. [HI the form of F NP (b) factorizes 
from the evolution. In both models there is no depen- 
dence of UPD on x at the initial scale Qq. As a results, 
we get the following set of initial moments: 



f^ s °\n,b, Q ) 

fr\n,b, Q ) 
fh vo \n,b,Q ) 



1 



71 + 1' 
1 

n + 1' 
0. 



(4.16) 



We remark that away from the chiral limit the separa- 
bility of the dependence on x and b no longer holds. In 
this case the initial conditions for the evolution are more 
complicated (they depend on b), but the analysis can be 
easily generalized to account for this case as well. 




FIG. 2: The evolution-generated UPD's for the pion for var- 
ious values of the transverse coordinate (from top to bottom 
b = 0, 1, 2, 3, 4, 5 and 10 GeV -1 ), plotted as functions of the 
Bjorken x. The evolution is made with the initial condition 
(ifcfll at Q = 313 MeV up to Q = 2 GeV. Solid lines - gluons, 
dashed lines - valence quarks, dotted lines - sea quarks. 



V. NUMERICAL RESULTS 

In Fig. [2] we present the results of our numerical cal- 
culation with the method using the Mellin transform. 
The initial conditions are for the pion in the chiral limit 
1)4. 16|) . holding at Qo = 313 MeV, and the evolution is 
carried up to the scale of Q = 2 GeV. The differential 
equations for the moments, H3.4I3.6J) . are solved numeri- 
cally for complex n— values along the Mellin contour, and 
subsequently the inverse Mellin transform l|3.9|l is carried 
out. The figure contains three families of curves, solid for 
the gluon, dashed for the valence quarks, and dotted for 
the sea quarks. In each family b assumes the values 0, 
1, 2, 3, 4, 5, and 10. Naturally, increasing 6 results in 
a decrease of the distribution, with the effect strongest 
at low x. At x close to 1 this effect disappears, which is 
explained in Sec. IVIIII We also note that at high values 
of b the distributions for the gluons become negative, re- 
flecting the change of sign of the Bessel function in the 
evolution kernel of Eq. 1)2. 4[l . As already discussed in 
Ref. [j"o]|. this poses no immediate physical problems, as 
the distributions in k± remain positive as primary ob- 
jects, and so are the physical cross sections. 

The results of Fig. are consistent with the findings of 
Ref. p^ |. where a different numerical method was used, 
as well as a different initial condition tested. 

In Fig. |31 we show the dependence of the evolution- 
generated UPD's on b at Q = 2 GeV and x = 0.1. 
The results are represented by squares for the non-singlet 
quarks, diamonds for the singlet quarks, and stars for the 
gluons. We note the much faster fall-off with b for the glu- 
ons than for the quarks, as expected from Eq. (|7.8|) . The 
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FIG. 3: Evolution-generated UPD's for the pion for Q — 
2 GeV and x = 0.1, plotted as functions of b. The evolution 
is made with the initial condition l|4.4[l at Qo = 313 MeV. 
The numerical results are represented by squares for the non- 
singlet quarks, diamond for the singlet quarks, and stars for 
the gluons, while the solid line shows the asymptotic formula 
17.51 for the case of non-singlet quarks. We note the much 
faster fall-of for the gluons than for the quarks, as expected 
from Eq. 17.81 . The quarks exhibit a long-range tail, accord- 
ing to the power- law formula 1 17.51 . As Q is increased or x 
decreased, the distributions in b become narrower, leading to 
spreading in k±. 



quarks exhibit a long-range tail, according to the power- 
law formula l|7.5|l . The solid line shows the asymptotic 
form for the case of non-singlet quarks from Eq. I|7.5I7.6|) 
(see Sec lVIIjl . which becomes accurate for b > 10 GeV -1 . 
As Q is increased further, or x decreased, the distribu- 
tions in b become narrower, leading to larger spreading 
in k±. For more details and plots concerning other nu- 
merical results see Ref. 
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FIG. 4: The low-fa expansion for the evolution-generated 
UPD's of the pion at b = 5 GeV" 1 and Q = 2 GeV 2 , plotted 
as a function of the Bjorken x. Solid lines - gluons, dashed 
lines - valence quarks, dotted lines - sea quarks. For each 
kind of parton the curves from top to bottom correspond to 
2, 4, 6, 8, 10, and 16 terms in the small 6-expansion. The 
initial condition for the evolution is provided by Eq. 14.41 . 



where B is the Euler beta function. We note that al- 
though at the level of the differential equation our ex- 
pansion is formally in Qb, the result l|6.2[l . together with 
the fact that r k is proportional to AQ fe CD (cf. (I A4|) ) . show 
that the expansion parameter is actually 6Aqcd- The 
rate of convergence of the method can be deduced from 
Fig. (@J. As we can see, the number of terms needed in- 
creases for increasing b and decreasing x. For x > 0.01 
eight terms in the expansion seems more than sufficient. 



VII. LARGE-6 EXPANSION 



VI. LOW-fe EXPANSION 



Since the anomalous dimensions (|B6IB7|) involve gen- 
eralized hypergeometric functions, they are not easy to 
use in numerical calculations. For that reason we con- 
sider the small-6 expansion, as well as the asymptotic 
forms at large bQ, presented in the next section. It is 
convenient to introduce the notation 



r k =r k (Q 2 ,Q 2 ) 



Q - dQ' 2 a(Q' 2 ) 2k 
Ql SvrQ' 2 



(6.1) 



The explicit form of functions is given in Appendix lAl 
Next, we apply Eq. I|C5H and find the following expansion 
in the non-singlet channel: 



f(n,b,Q) 
f(n,b, Qo) 

x [B(2fc,n- 



(0) 

e T„,jvs r » 



exp 



2^fe / ^l-fe 



fc=l 



(-b 2 ) fc 4 



B(2fc,n + 3)] r k ] 



(6.2) 



Appendix [D] contains asymptotic forms of the gener- 
alized hypergeometric functions and of the 6-dependent 
anomalous dimensions. These expressions hold in the 
limit of large Qb at n kept fixed. The formulas are of 
great practical importance in the present study, since 
they are much simpler to implement in numerical cal- 
culations than the generalized hypergeometric functions 
appearing in Eq. HB6IB7jl . Actually, Eq. (|D3ID6|I should 
be used whenever Qb is larger than about 10|n|. Since 
in practical problems the scale Q may be as large as the 
mass of the IT-boson, there is a frequent need to use the 
asymptotic expressions (|D3ID6|) . Also, if the UPD's in 
the transverse momentum are needed, one has to carry 
back the Fourier transform from the &-space to the k±- 
space, which involves all values of b. 

We start with the non-singlet case of Eq. I|3.4|l . With 
the help of Eq. I|D3(I , where at the leading order in 1/ (Qb) 
we drop the oscillatory parts, we may write l|3.5|) at large 
Qb: 



/ns°VAQ) 
fWi\n,b, Qo) 



oVo+Vin 



(7.1) 
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where 



and 



yi 



Vo 



Z ^QCD 



_4C F ( _^-r_ 3/2 (Q3> Q 2 ) + — log 



, foA QCD , 3 
log —^-+7 — 



l , Q 

p { 3 > Qo 

MQlQ 2 ) 

for Q 2 < n 2 c , (7.2) 



and 



1)1 



-4C F (^%Hr_3 /2 (Q^^) 



^-3 /2 (^,o 2 )-^_3 /2 (g^Q 2 



yo 



-4C f 



2A 



QCD 



2A 



1 



QCD ,„ 2 2\ 

3" 



6A 



log 

z 

log 

S 2 



r (Mc,Q ) + 



7 "4 



for Q 2 > n\. (7.3) 



Next, we take our initial condition 1(4.16(1 and use the 
inverse Mcllin transform 



dn 
2Td" 



n + 1 



: " - yi Q [yi +log- ) ■ (7.4) 



The condition provided by the theta function means that 
the formula can be used for x < exp(yi). For negative y\ 
this means that the validity is limited for x not too close 
to 1. This, however, has been already tacitly assumed, 
since the asymptotic expansion holds for fixed values of 
n, hence cannot describe x in the vicinity of 1. The above 
formulas lead to the following large-6 form of /n™ 1 m the 
x space: 



fSt{x,b,Q) = x 
3~ 



fbA 



QCD 



x exp 



27- 



1 



10 
i2 ,n2 



V 2 

ro(Ql,Q 2 ) + ^r_ 3/2 (Qi(fi 

for Q 2 < (x 2 c , (7.5) 



A, 



x A" 



-8C F r ( M ;?,Q 2 ) / Mc 



-4CV/AT 



-BC F r {Ql,p. 2 c ) 
QCD 



x exp 



27 



rn(QlQ 2 ) + ^r_ 3/2 (QlQ 2 ) 



for Q 2 > /i 2 , 



(7.6) 



We note a few facts: the non-singlet UPD of the 
pion is for large Qb linear in x for x not too close 
to 1, with the slope decreasing with b as a power 
law (we can neglect here the small correction due to 
the last term in Eq. 1(7.517.6(1 . The exponent of b is 
— 8Cfto(Qo, Q 2 )- The overall constant is also deter- 
mined. Note that linearity with x at not-too-large x is 
seen in Fig. [21 for the valence quarks (dashed lines) at 
b = 10 GeV- 1 . Numerically, at Q = 2 GeV we find 
for low x that f NS (x, 6 = 5 GeV" 1 , Q = 2 GeV) = 0.73x 
and / NS (a;, b = 10 GeV" 1 , Q = 2 GeV) = 0.35a;, in ac- 
cordance with the exact calculation of Fig. [3 

The asymptotic form ((7.517.6(1 works very efficiently 
in practice. For the case of the non-singlet quarks this 
can be seen from Fig. [3J where for Q — 2 GeV and 
b > 5 GeV -1 there is virtually no difference between the 
exact numerical calculation and the asymptotic formula. 

The power-law behavior in b shows that the evolution 
generates a rather weak behavior at large b. Numerically, 
at Q = 2, 4, and 100 GeV, the power of b is, respectively, 
-1.12, -1.29, and -1.75. This means that the large-6 be- 
havior for the non-singlet quarks is controlled by the ini- 
tial profile F NP (b), which in chiral quark models of the 
pion has an exponential fall-off, rather than by the QCD 
evolution. 

Now we pass to the discussion of the singlet case of 
Eq. 1)3.6(1 . In the large-Qfe limit the leading part of the 
matrix T n of Eq. ((3.7(1 becomes 



T„(Q6) 



4C F log^ 

4iV c log^ 



(7.7) 



From this form, using methods as for the non-singlet case 
above, we infer that the dependence on b is asymptoti- 
cally 



f s {x,b,Q) 
fa(x,b, Q) 



b -SC F r (Ql,Q 2 )^ 
b -SN c r (Q 2 ,Q 2 ) _ 



(7. 



Thus the singlet quarks fall off at the same rate as the 
non-singlet quarks of Eq. ((7.517.6(1 , while the gluons drop 
significantly faster, as C/ = 4/3 and N c = 3. This be- 
havior is clearly seen in Fig. [3J We note that due to 
the complications of Eq. 1(3.7(1 a more detailed analysis 
yielding pre-factors, such as the one for the non-singlet 
case presented above, is more difficult in the present case, 
hence we do not pursue it further here. 



We end this section with a couple of remarks concern- 
ing the observed long-range nature of the tails in b. Our 
initial non-perturbative profiles F NP (6) drop exponen- 
tially, therefore suppress the tails generated by the evo- 
lution. This means that the large-6, or low-fcj_ behav- 
ior is controlled by non-perturbative effects. The larger 
negative power in fa, Eq. I|7.8|) . explains the the faster 
shrinkage of gluon distributions in the &-space, or faster 
spreading in the fcj_-space. 



VIII. BEHAVIOR OF / NS AT x 




According to standard properties of the Mellin trans- 
form, the limit x — > 1 is obtained form the large-n be- 
haviour of the anomalous dimensions. Using the asymp- 
totic form B(n,m) — ► r ^ in Eq. I|C4|I . or the explicit 
form of the anomalous dimensions (|B6IB7|) . one obtains 
that 



Jn,Ns(Qb) - 7 Svs = 2C. 



b 2 Q 2 



+ 0(l/n 3 ). (8.1) 



We also need the large-n expansion of 7 



(o) 

n,NS> 



which is 



ln,NS = (logn + 7 - 3/4 + R{n)) , (8.2) 

where R(n) — J2kLi Ckn ~ k ■ Thus, according to Eq. H3.5|) . 
we have the asymptotic form 



= n~ 8CFVo( - 



fSs\n,b,Q 2 

x e -8CFro(7-3/4+/?(n)) e 2CF6 2 ri/n 2 +C>(l/n 3 ) 



(8.3) 



Using the initial moments (|4.16l) and expanding the ex- 
ponentials in Eq. (|8.3|l we obtain 



f^ s °\n,b,Q 2 ) 



-SC F r 



V fe=l / 



x {l + 2C F b 2 r 1 /n 2 +C(l/n 3 )) 



(8.4) 



Next, we use the formula for the inverse Mellin transform 

- A r(A,-wio g ^ 



dnx n 

c n + w 



1 



T(A) 



1 \A 



AT (A) 



(8.5) 



which after a straightforward algebra leads to the equa- 
tion 



2C F b 2 r 1 {l-x) 2 
1_ (1 + 8C F r )(2 + 8C F r ) 
+ 0{{l~xf). (8.6) 



The terms with coefficients c' k do not enter at the leading 
order in (1 — x). The behavior of Eq. I|8.6|) agrees with the 



FIG. 5: Model prediction for the integrated valence quark 
distribution of the pion, evolved to from the initial condition 
14.1611 to the scale of Q — 4 GeV, confronted with the E615 
Drell-Yan data HI. The behavior at x -» 1 is (1 - x) 1 - 29 . 



behaviour of Fig. [3 where close to x — 1 the departure 
of the curves with finite b from the curve with 6 = 
becomes very slow, as it is suppressed by (1 — x) 2 . 

We also obtain that at x — ► 1 the integrated non-singlet 
distribution behaves as 



/ns(z,0,Q 2 ) 



e 2C F (3-4 7 )r 

r(l + 8C F r ) 



(1 



\SC F r 



•7) 



This agrees with the fact that a function which originally 
behaves at x — * 1 as fisis(x, 0, Qo) — ► (1 — x) p evolves into 

M 



frs(x,0,Q 2 )^ (l~x) p - 



00 



■ log : 



(Qo) 



In our approach the integrated function at Q 
p = 0. Numerically, we find that 



/ns(M,(2 GeV) 2 ) 
/ns(x,0,(4 GeV) 2 ) 



1.15(1 -a;) 1 - 13 , 
1.08(1 - x) x - 2Q . 



Qo has 



(8.9) 



Note that although with the DGLAP evolution the 
Brodsky-Lepage counting rules for the behavior at x — > 1 
are clearly disobeyed, our numerical predictions agree 
within experimental uncertainties with the experimen- 
tal data, including the region very close to x = 1. Fig.0 
confronts our results evolved to the scale of Q = 4 GeV 
to the E615 experimental Drell-Yan data [58| which cover 
the large- a; region. In view of the simplicity of the present 
model, the quality of this comparison is impressive. 



IX. BEHAVIOR AT x 







The low- a; behaviour of the inverse Mellin transform 
is encoded in the closest singularities to the integration 
line. We start with the non-singlet case of Eq. (|3.4J) . 
Using the pole-residue expansion of Appendix lEl we find 
that the the closest singularity is at n = — 1, with the 
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residue — 4Cf Jo(Qb). With the help of the expression 
for the inverse Mellin transform, 



dn x~ n _ I x/ (2^/alogi), a >o 

2m n + 1 " I „ T _ro L i ' 



we find that 



where 



.4 



a;Jo(2,/alogi), a<0 



(9.1) 



xI (2yJC F Alog±), A>0 
xJ (2Jc F Alog±), A<0 



dQ 2 
2irQ' 



■a{Q 2 )J {Qb). 



(9.2) 



(9.3) 



For the singlet case we retain the closest singularity at 
n = and rewrite Eq. (|3.7|) in the form 



where in the last line we have used the asymptotic form 
of the Bessel functions. With this result we find that the 
unintegrated gluon distribution at x — > behaves as 



/ G (n,6,Q)~exp^2iV c AlogiJ , A > 0, (9.7) 

with A provided in Eq. (|9.3|) . If a < in Eq. 19.6[> . then 
the Ik functions above are replaced with the Jk func- 
tions, and the asymptotics changes the character from 
exponential to oscillatory: 



f G {n,b, Q) 



1 loe 



-.4 



cos j .1 log - 



1 - 



^) S in(2 V -Alogi 



A < 0. 



-Atagi(l-^ 



(9.8) 



exp 



fs{n,b,Q) 
f G {n,b,Q) 

Q2 dQ' 2 a(Q' 2 )J (Q'b) 



(9.4) 



ttQ' 2 





C F A. 



/s(n, 6, Qo) 
fG(n,b, Q ) 



where we have used Eq. |E6J), and could drop the V sym- 
bol since in the present approximation the matrix in the 
exponent does not depend on Q' . After some straight- 
forward algebra we obtain 



fs(n,b,Q) 
f G (n, b, Q) 



fs(n,b,Q ) 



/sM,Qo)-^(e— -l) 



+ fG(n,b,Q )e~ 



N c 

2N C A 



(9.5) 



For b = the result (|9.7|l is consistent with the 
double-leading-logarithmic approximation (DLLA) for 
the DGLAP equation where one obtains, with con- 
stant a, 



xg(x, Q) ~ exp 



N, 



Q 2 



l 



— a log log - 



(9.9) 



See, e.g., the review in Ref. |6(j. Our formulas i|9.2|) and 
(19.7(1 are generalizations of this behavior for the uninte- 
grated distributions evolved with Eq. 1)2. 4f> . 

The pole-residue expansion of Appendix lEl is good for 
not-too-large Qb. This limitation, at any fixed x, carries 
over to Eq. I|9.2|) Numerically, we find that at x = 0.1 
Eq. (|9.2|) is valid for Qb < 5. For higher values corrections 
from further residues should be included. 



The equation for fs shows the inadequacy of retaining 
for this case the singularity at n = only, as in the 
considered limit of x — > the singlet quarks are controlled 
by the singularity at n = —1. For the case of gluons 
we use the initial condition H4.16|) . We need the inverse 
Mellin transform 



dnx " — 

c n 

oo 

k 



1 



fc=0 



— exp 
1 \n 

fe/2 



fog h 



Z fc (2Walog-) 



exp (2Walogi 



47r A /alogi 1 



a>0, (9.6) 



X. EVOLUTION OF (k 2 ± ) 

The average transverse momentum squared is a con- 
venient measure of the width of the UPD's. Due to the 
factorization of the initial profile, Eq. 1)2. 6fl . (k 2 ) decom- 
poses into two terms: the contribution from the initial 
profile F NP , gives the width at the initial scale Q = Qq, 
and the piece (fc|,) cvo i, entirely to the evolution and in- 
dependent of the profile F NP , 



(k 2 ± ) = (fci) NP + (feDevol, 



(10.1) 



(ki) 



NP 



(fc_|_)evol 



= -4 



= -4 



dF NP {b)/db 2 



b=0 



F NP (b) 
df cvo \b lX ,Q)/db 2 



f ev ° l (b,x,Q) 



6=0 
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The contribution (fcf^NP has already been discussed in 
Sec. I1VI hence here we analyze the term generated by 
the evolution. 
Let us denote 



m df° vo Hb,x,Q) 



(10.2) 



6=0 



In the Mellin space, the equations obtained by expanding 
Eq. (|3.4|) up to first order in b 2 around b = read 



dQ 2 



+ Q 2 7i']vs/Ns(«,0,Q)) , (10.3) 

and 

djjp (n, Q) _ as(Q 2 ) ( (o) Ai) ( Q) 

+ l^ qG f i G ) (n,Q)+ (10.4) 
Q 2 [l ( nl q h(n, 0, Q) + J^ G f G (n, 0, Q)]) , 



8nQ 

QVnhM^ °- Q) + Inhcfain, 0, Q)] 



which form a set of ordinary inhomogeneous differential 
equations. Since at the scale Qo all the width is by con- 
struction generated by the initial profile F, the initial 
conditions for Eq. H10.3I10.6JI are 



/, (1) (n,Qo) = 0. 



>..... e (10.6) 

For the non-singlet case we have the solution 



(!) 



(10.7) 



In the singlet channel we carry the analysis numerically. 

Next, we pass to the x-space via the numerical in- 
verse Mellin transform. The results for the dynamically- 
generated root mean squared radius of the pion are shown 
in Fig. for various values of x. In confirmation of the 
results of Ref. [llj , we note that the k± width increases 
with Q for all parton species. The width for the glu- 
ons (solid lines) is larger than the width of the non- 
singlet (valence) quarks (dashed lines), and the singlet 
quarks (dotted lines). With the log-log scales of Fig. 
the slopes of the plotted lines become to a good approx- 
imation equal to one another at large Q 2 . 

With the help of previously-derived expressions for the 
behaviour of /ns near x = and x = 1 we may obtain the 
following expressions (fej^NS near the end-points. From 



Eq. H8.8f> we have at x — > 

hy-AC F r log a;) 
ib (\Z-4GVo log a) 



lh 2 \ evo1 



Cf log x 

ro 



-n 



Cf log x 

To 



-r\. 



(10.8) 



At large Q 2 the leading behavior is 



(k 



2 \evol 
_L/NS 



2f3^C F log^a(Q 2 ) ( 

\ w _M) 8?r 

N g a (Q 2 ) 



(10.9) 



i.e. up to the log log Q 2 corrections the spreading pro- 
ceeds with a(Q 2 )Q 2 . At x — > 1 we find from Eq. IjS.fijl 
that 



2 ]cvo i 2GV(1 -x) 2 n 

[ ±/NS (l + 8C F r )(2 + 8C F r ) 



(10.10) 



At large Q 2 the leading behavior is 



(k 



2 \evol 
NS 



^ 4)2 (1 



HQ 2 ) 



64ttC f 



log -Ml 



aQ 2 



Q 2 . (10.11) 



Again, the growth is, up to the log log Q 2 corrections, 
proportional to a(Q 2 )Q 2 . 

For the gluons the same asymptotic behavior of 
(^i)c o1 follows from Eq. (|9.6(l . Thus, to summarize, all 
UPD's grow at large Q as Q 2 a(Q 2 ), in accordance to the 
behavior in Fig. 

Interestingly, it can be noticed from Fig.|Hlthat at Q — > 
Qq the k± width for the gluons does not vanish. In 

this limit both /§ vol (a;, 0, Q) and /q\x,Q) vanish, as is 
obvious from Eq. (|4.2I10.6|I . Thus one has a 0/0 limit. 
From Eq. H2.4I10.3|) with the initial condition I|4.3I10.6|) 

one can easily obtain that 



km (k 2 x yj° l = Qo 

Q— »vo 



J x dzP Gq (z) 



(l-zf 



Gx d 



S x dzP Gq {z)\ 
21x 2 - 18x log x- 10a; 



(10.12) 
G 



3(x 2 — 2x log x + x — 2) 

which is positive for x E [0, 1) and equal for x = 1. On 
the other hand, since for the quarks 0, Q) ^ 

(^±}ns\s vanish at Qo- 

In phenomcnological applications it is sometimes use- 
ful to have a simple formula characterizing the discussed 
behavior. In the range 2 GeV 2 < Q 2 < 10000 GeV 2 and 
0.005 < x < 0.8 the following simple-minded interpolat- 
ing formula works to within a few per cent: 



((kl 



2 \ovol 



1/2 



A<(log-) 

X 



Q 2 



0.35+0.004 log ■ 



A 2 
QCD 



(10.13) 
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FIG. 6: The rms transverse momenta of UPD's of the pion for 
x — 0.01, 0.1, and 0.5, plotted as functions of the renormaliza- 
tion scale Q 2 . Solid lines - gluons, dashed lines - non-singlet 
quarks, dotted lines - singlet quarks. 



where i =NS, S, or G, and 



A G {y) 



-0.017y 1/2 
-0.021y 1/2 
-0.016?/ /2 



0.1132/ - 0.057y 3/2 + O.OlOy 2 , 
0.120y - 0.059?/ 3/2 -I- 0.009j/ 2 , 
0.150y- 0.075?; 3/2 + 0.011y 2 . 

(10.14) 



The power of Q 2 of 0.35 in Eq. (|10.13|l . rather than 1/2 
corresponding to (k\)^ vo1 ~ Q 2 a(Q 2 ), compensates, in 
the chosen range for Q, for the logarithmic corrections. 
We note that Eq. I|10.13|) holds for the pion with the ini- 
tial conditions 1)4. 16|) provided by the chiral quark mod- 
els. 



XI. FORMAL LIMITS FOR OTHER INITIAL 
CONDITIONS 



becomes 

dn 
(j 2iri 



pVo+Vin / | 

n + a \ ' 



and expressions (|7.5I7.6[1 are modified accordingly, with 
x replaced by x a , and y\ in the exponent multiplied by a. 
The formulas (|7.8|l remain valid. Therefore the power-law 
behavior at large b is independent of the initial condition, 
and is &-8CWQ^,Q 2 ) for the quar k s an d ^sN^oiQlQ 2 ) 

for the gluons. 

In the limit of x — > 1 the only difference is the appear- 
ance of the extra power of (3 in the in Eq. I|8.7|l . Thus, 
the UPD's for finite b approach the b — case of the 
integrated distributions as (1 — x) 2 . 

In the limit of x — > we need generalizations of the 
Mellin transforms of Eq. (|9.1I9.6[) for a^l. These are 



dnx ™ 

c n + a 



exp 



= xj2(-i) k (a-iy 



fc=0 



log; 



(11.2) 



fc/2 



xexp 



(2^/alogI 



a - l)W47rWalogi 



7 fe (2\/alog-) 



a^l, a > 0, 



and 



dnx n exp ( — 

, • n + a \n 

fc/2 



(11.3) 



= Ec-d* 



fc f lo s^ 



fc=0 



7 fe (2 A /alog-) 



exp (2 Jalogi 



ai/47r^/alo, 



a ^ 0, a > 0. 



The analogs of Eq. (|9.2l9.7ll0.8ll0.9jl follow straightfor- 
wardly. In particular, the Q 2 a(Q 2 ) large-Q behavior for 
all parton distributions is preserved. 



XII. CONCLUSIONS 



Certain formal results listed in this paper, such as the 
formulas for the non-singlet quarks, (|7.51 17.61 18.71 18.81 
19.21 110. 8110. 9fl . are specific to the evolution with the ini- 
tial condition following from the chiral quark models, 
Eq. I|4.4I I4.l6|) . However, these results can be easily gen- 
eralized. Note that most of the popular parameteriza- 
tions of initial conditions, such as those of Refs. |52u53|. 
involve linear combinations of x a (l — x)@ . It is under- 
stood that the factorization in the initial condition be- 
tween x and b variable holds, as assumed in Ref. [T^ . 

For the case of large-6 asymptotics, the relevant for- 
mula is (|7.4|) . With the initial condition x a (l - x) p it 



We have presented a new method of solving the 
Kwiecihski equations for the leading-order QCD evolu- 
tion of unintegrated parton distributions. The method is 
based on the Mellin transform and parallels the standard 
analysis of the DGLAP equations. Our main results are 
as follows: 

1. We have found analytic forms of the 6-dependcnt 
anomalous dimensions, expressed through hyperge- 
ometric functions, which allowed us to study formal 
aspects of the equations and their solutions, e.g. 
the asymptotic forms of the evolution-generated 
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UPD's at large 6, or at x — > and x — ► 1. We 
have also demonstrated that the proposed numeri- 
cal method is fast and stable. 

2. The numerical work can be simplified if low-6 or 
large-6 expansions are used. 



EURIDICE grant number HPRN-CT-2003-00311 is ac- 
knowledged. Partial support from the Spanish Ministc- 
rio de Asuntos Exteriores and the Polish State Commit- 
tee for Scientific Research, grant number 07/2001-2002 is 
also gratefully acknowledged. 



3. At large b the evolution-generated 6-dependent 
UPD's exhibit power-law fall-off, with the magni- 
tude of the exponents growing with the probing 
scale Q, cf. Eq. (|7.5I7.8|I . The fall-off is steeper 
for the gluons than for the quarks. 

4. At x — > we have found generalizations of the 
DLLA behavior, cf. Eq. jQQ ) We have also 
shown that for large b the solution for the valence 
UPD of the pion grows linearly with x for not too 
large x, and the slope decreases with b as a power 
law. 

5. At x — > 1 the evolution-generated 6-dependent 
UPD's approach the integrated distributions as 
(1-x) 2 . 

6. Our numerical results fully confirm the finding of 
Ref. |l6j |. where a different numerical method was 
used. We find the spreading of the k± distributions 
with the probing scale Q, with the effect strongest 
for gluons and increasing with decreasing x. We 
have also shown that the widths (k 2 )^ 01 in all 
channels i increase at large Q 2 as Q 2 a(Q 2 ). 

7. For practical purposes in possible phenomenolog- 
ical applications, we have parameterized (fc^)™ oi 
with a simple formula which works with accuracy 
of a few percent. 

Although the specific study of this paper was devoted 
to the pion with the initial condition following from the 
chiral models, and several of the more detailed analytic 
formulas were specific to this case, the developed method 
is general and can be applied to any initial form of the 
data. In particular, it can be used with the GRS [H^ 
or GRV |53j | parameterization supplied by a profile in b, 
such as already studied in Ref. [Tg • The formal results of 
Sec. IXII are general for a wide class of initial conditions, 
suitable for both the pion and the nucleon. 
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APPENDIX A: ELEMENTS OF THE 
PERTURB ATIVE QCD EVOLUTION 



up to the scale fi 2 — 4 GeV 2 and four active flavors above. 



We use the LO QCD evolution with three active flavors 
Therefore a = g 2 /{Air) is given by 



4tt 



a(Q 2 ) 



ti 3) log 
4tt 



JJ1 



l QCD 



A, 



A, 



QCD 



with ^ Nf) = 11 - 2AT//3 for N c 
N c denote the number of flavors and colors, respectively. 
Along this paper we take 



Q 2 < & 
Q 2 > /& 

(Al) 

3, where Nt and 



Aqcd = 226 MeV, 

as was done in Refs. pil I4 H lil Ul\ 

scale A4 ensures matching at Q 2 - 



(A2) 



= 9, (3^> = 25/3, and A 4 = 189 MeV. 
The functions , defined in Eq. (jA"|l , have the explicit 
form 



-Al) 



The value of the 
Numerically, 



ro(QlQ 2 



1 



m 
1 



(3) 



log 



■\og(Q 2 /A 2 QCD Y 
log(Q2/ A 2 QCD ) 



ro(QlQ 2 ) 



ro(Ql,H 2 c ) 



\ 2k 



log 



«(Q§) 
a(Q 2 Y 

1 



1 



(4) 



(4) 



log 



log 



Q 2 < Mc, 



log(Q 2 /Al) 
log(^/A|) 



a{Q 2 Y 

Q 2 > Mc 



(A3) 



and for k 7^ 
r k (QlQ 2 ) = 



r k (QlQ 2 ) - r k (Q 2 ,£) + 




A 



2/,- 



(4) 



Li 



Af 



" 1 % 



(A4) 



Q 2 > A- 
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Above we have used the indefinite integral 

U 3 v 



-A-uf^.*-!,*,. 



<2 2 log(<3 2 /A 2 ) 

where the logarithmic integral is 

Li(x) = / dt/logt. 
Jo 

At large Q 2 



(A5) 



(A6) 



fclog^r fc 2 log 2 



(A7) 



The functions P a t,(z) are the LO splitting functions 
corresponding to real emission, i.e.: 



P qG (z) = Nf [z 2 + (l-z) 2 } 

PoM = c,l±i^!, 

P GG {z) = 2N C 



(A8) 



2 1 — z ,„ . 

+ + z(l-«) 

1 — z z 



with C F = {Nl - \)/{2N c ) = 4/3. 



APPENDIX B: 6-DEPENDENT ANOMALOUS 
DIMENSIONS 



We introduce 



Q 2 b 



2 h 2 



(Bl) 



as well as the anomalous dimensions for the b — case 
of the DGLAP equations, where for the non-singlet we 



have 



ln,NS 



2C f 



1 + n 



4H n ,(B2) 



while for the singlet 



(B3) 



7& = -4C P (-- 2 



1 + n 2 + n 3 + n 
1 



<y (0) - 2A^ (-3-- 



n 1 + n 2 + n J 
4 



n 1+n 2+n 3+n 
+ 4JT n ). 

The symbol i? n denotes the harmonic number 



B fc!* r(n + i) +7 ' 



(B4) 



which is a meromorphic function in the complex n 
variable, with poles located at negative integers n = 
— 1, —2, —3, ... and residues equal to — 1. 

Below we list the anomalous dimensions for the mo- 
ments of the unintegrated parton distributions in the b 
space, defined in Eq. I|3.2|l . The formulas follow from the 
basic analytic integral 



r(2 + n + v) 



r(i + M )r(i + I /) Jo 

. 1 + u 2 + fi 2 + fj, + v 3 + + v 



dyy^(l-yyj (2^1y)) 
-u) 



(B5) 



and relations among the generalized hypergeometric 
functions. For the non-singlet case we have 



ln,N S (Qb) = 7$ S + (1+ — 2 + n) 



4C F 



-3 - 2n + 2 (2 + n) X F 2 {1 *±±, *±^; -u) 



- iF 2 (-; -u) + 2u 3 F 4 (l, 1, -; 2, 2, -u) 



(B6) 



whereas for the singlet case 



ln,qq{Qb) = -fn,Ns{Qb), 
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ln,qG{Qb) 



- J°) 



4N 



+ (l + n)(2 + n)(3 + n)(4 + n)(5 + n) [ ~ ((4 + ^ (5 + " } <" 4 ~ " (3 + "> 
+ (2 + „) (3 + „) 2 F l{ l i + -u) - 2 (3 + n) xF a (|; i±^; -«) 



/1E ,,3 4 + n5 + n , 
4lFa <2 ; — — ;_tt) 



„,56 + n7 + n 
24hiF 2 (-;— ; -u) 



ln, Gq {Qb) - 7i°G g + n (1 + n) (2 + + n) (4 + n) 



4C F 



( (I i nU2 n) (3 + n) (4 + n) ^(i; ±±^, - U ) 

+ ! 3 + , ) ( I i „ 1 ( 4 + n (3 + „) - 2 i±^; -u)) + 12 u ^ -«) 



1 



1 



fKi; l + § , | + f ; -«) iF 2 (|; 1 + 1, | + 1; -«) 



n 1 + n 2 + n 3 + n 
F 2 (|;| + f,2 + f;- U ) 2 1 ^ 2 (|;2 ' " ' : "■ 



1 +n 



n + 71/ 



2 ' 2 1 2 ' 



-«) 12iiiF 2 (|;3 + 5,| + §;-u) 



(1 + n) (2 + n) 



(1 + n) (2 + n) (3 + n) (1 + n) (2 + n) (3 + n) (4 + n) (5 + n) 



^3^4(1, 1,|; 2, 2, § + H,2+f;-u) 



(1 + ti) (2 + n) 



(B7) 



One may verify that the analyticity properties in n of we arrive at the expansion formulas 
the anomalous dimensions (|B6IB7|) are the same as for 
the 6 = case of (|B2IB4ll . 



Jn,Ns(Qb) = ln,qq{Qb) = J$ qq 



APPENDIX C: EXPANSION OF ANOMALOUS 
DIMENSIONS AT LOW bQ 

We may expand in the anomalous dimensions in powers 
of Q 2 b 2 , 



which yields 



(CI) 



ln,NS "n,qq 



ln,qG 



2C F (n 2 +5n + 7) 
~ (ti + 1) (ti + 2) (?i + 3) (ti + 4) : 
2N f (n 2 + 3n + 14) 



ln,Gq 



in.GG 



(n + l)(n + 2)(n + 3)(n + 4)' 

2C F (n 2 + 7n + 24) 
(n + l)(ra + 2)(n + 3)(n + 4)' 

+ 5)(n 2 + 5n+ 16) + 120] 
7i (n + l)n + 2)(ti + 3)(n + 4)(n + 5) 



(C2) 



More generally, introducing the Euler Beta function, 
B(x,y) = T(x)T(y)/T(x + y), and applying the series 
expansion of the Bessel function, 



j (x) = £ 



2\k 



k=0 



(~* 2 ) 

2 k k\ 2 ' 



(C3) 



k=l 



2l.2\kAl-k 



(-Q 2 b 2 )H 



ln,q G (Qb)=l { Z G -N f J2 



[B(2/c,n + l) + B(2fc,7i + 3)] 
{-Q 2 b 2 ) k 4 1 - k 



k=l 



k\ 2 



x [B(2fc + l,7i+ 1) - 2B(2fc + l,n + 2) 
+2B(2fc + 1,71 + 3)] 



7n, Gg (Qb)=7i° G9 -^E 



(-Q 2 b 2 ) k 4 1 - 1 



k=l 



x [2B(2/c + 1, ri) - 2B(2/c + 1, n + 1) 
+B(2fc + l,n + 2)] (C4) 



7n, GG (g&) = 7i 0) 



GG 



2u2\kAl-k 



k=l 



(-Q 2 b 2 )H 
W 



x [B(2fe, ri + 2) + B(2fc + 2,n) + B(2fc + 2, ri + 2)] 



APPENDIX D: ASYMPTOTICS OF THE 
ANOMALOUS DIMENSIONS AT LARGE bQ 

We may use the asymptotic forms of the generalized 
hypergeometric functions appearing in Eq. (|B6IB7|) . One 
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has 165 



iF 2 (ai;b 1 ,b 2 ; 



Q 2 b 



2u2 



r(6i)r(6a) 



T{h - ai)r(6a - a x ) VQ 2 ** 2 



r(bi)r(6 2 ) 



, 4iV> / Q 2 6 2 11 

7GG(n, Q6) = -g^ + 4N C (log + 2 7 - — 

2n + 2 4T(n)(Q6)-"-3 ,2n + l 
+ coaf-^-Q*) 



{cos(Q6 — TTCl) 

+ 5 sin(Q6 ~ 7rCl) }(ot 



[r(n)-20r(n+l)](Qb)-"-2 

77* 77 \ 

x [cos(-tt - Qb) + sin(-7r - Qb)]) + 



(D6) 



ci = - (h + b 2 - a\ - i 



1 



(Dl) 



c 2 = - (12a 2 - 8(6i +b 2 + l)oi - 4(6i - 6 2 ) 2 



+8(6i+6 2 )-3) 



and l6£ 



3 n + 3n + 4 Q 2 6 2 

3^4 1, 1, 7T, 2, 2, 



2' ' ' 2 
4(n + l)(n + 2) 



2 ' 



nQ 3 b 3 
8r(n + 3) 



Qb[nlogQ- mp°(n) — 1] + n' 



2?r 



cos ( Q6- 2n + 3 ?r ) + 



(D2) 



Then the following asymptotic expansions for the in- 
dependent anomalous dimensions hold: 



7ns (n, Q6) = 7 OT (n, Q6) = 4C F (fog + 2 7 
2n + 2 r(ri+l)(Q6)-"-3 



Qb 4%/2^ 
246Q cos( 2 " 4 '' 3 7r- Qb) 

2tj 4- 3 

-(12n + 13) sin(^— 7T - Q6) 



(D3) 



7 9 G(n,g6) =AN f 



1 r(n + l)(Q6)-™-2 



Q6 



8^ 



(-12n- 



71 7T 

-ll)cos( T -Q6) 



-(12n + 8Qb + 11) sin(— - Qb) 



(D4) 



The ellipses denote terms subleading in ^g. The above 
formulas assume that n is kept fixed. In actual applica- 
tions, such as numerical programming of the generalized 
hypergeometric functions, it is practical to switch from 
the general formulas (|B6IB7|) to the asymptotic expres- 
sions (|D3ID6jl when Qb > I0\n\. 



APPENDIX E: POLE-RESIDUE EXPANSION OF 
THE ANOMALOUS DIMENSIONS 

For 6^0 the analytic structure of the b— dependent 
anomalous dimensions remains the same as for 6 = 0. 
This can be seen by expanding the Bessel function in the 
integrand of Eq. i|3.2|) as a power series around 2 = 0, 
which yields 



ln,ab{Qb) 



(El) 



r 1 °° 1 
4 / dz^^[4 k \Qb)(-Qb) k z n+k -l 



Pab(z). 



k=0 

Applying the trick 



oo 1 

1 = J (Q) = -4j2^4 k) (Qb)(-Qb) 



k=0 



(E2) 



we find the expansion involving index-shifted anomalous 
dimensions at 6 = 0, namely 



fe=0 



{-Qbf tW/ 



ln,ab{Qb) = ^^^J W (Q6) 7 „+fc,a6(0). 



(E3) 



Using the explicit expression for the anomalous dimen- 
sion this series may be rewritten as a pole-residue expan- 
sion 



7«,«6(Q6) = E-jTir- 



(E4) 



k=0 



(4n + 8Q6- l)cos(— - Qb) 



-(An- 



71 IT 

l)sin(—-Qb) 



(D5) 



In practice this means that the Mellin contour used in 
the case of b = can be used in the 6^0 case as well. 
In the non-singlet case the first few residues read 



R™(A) 

rT(a) 



-4CWo(A), 
AC F (J (A) 
-2C f {-(A 2 - 



(E5) 



-AJ^A)), 

4)J (A) + 3AJi(A)), 
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while in the singlet channel R\ q = Rf s and 

Rf(A) = -4N F J (A), 

Rf(A) = -4N F (-2J (A) + AJi(A)), 

R Gq (A) = -8C F J (A), 

R Gq {A) = -8C F (-J (A)+AJ 1 (A)), 

R GG (A) = -8N C J (A), 

R GG (A) = ~8N C (-J (A)+AJ 1 (A)). (E6) 



The pole-residue expansion controls the behavior of the 
solutions of Eq. 1)2.4(1 at low x. Since the subsequent 
residues carry powers of A n — (Qb) n , the expansion can- 
not be used for Qb too large. 
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